R1 Indexes Cortical Myeloarchitecture
Read in Data for Myeloarchitecture Characterization
Final study sample
participants <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/sample_info/7T_MP2RAGE_finalsample_demographics.csv") #215 scans from 140 subjects. csv generated by /sample_construction/finalsample_7Tmyelin.RmdGlasser regions list
glasser.regions <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/HCPMMP_glasseratlas/glasser360_regionlist.csv")
glasser.snr.exclude <- read.csv("/Volumes/Hera/Projects/corticalmyelin_development/Maps/SNR/glasser_SNR_exclusion.csv")Depth-dependent R1 in HCP-MMP regions
#R1 by region matrices at 7 cortical depths
myelin.glasser.7T <- readRDS("/Volumes/Hera/Projects/corticalmyelin_development/BIDS/derivatives/surface_metrics/depthR1_glasseratlas_finalsample.RDS") #generated by /surface_metrics/surface_measures/extract_depthdependent_R1.Rordered_depths <- c("depth_7", "depth_6", "depth_5", "depth_4", "depth_3", "depth_2", "depth_1")
depth_colorbar <- c("#004A38FF", "#006458FF", "#0093B0FF", "#8AABD5FF", "#C0BFE3FF", "#E0D5ECFF", "#F2E7F4FF")Regional average R1 in HCP-MMP regions
#Calculate across-depths average regional R1
regionalaverage.myelin <- do.call(rbind, myelin.glasser.7T) #R1 at all depths
regionalaverage.myelin <- regionalaverage.myelin %>% select(contains("ROI"))
regionalaverage.myelin <- colMeans(regionalaverage.myelin) %>% as.data.frame() %>% set_names("mean.R1")
regionalaverage.myelin <- regionalaverage.myelin %>% mutate(orig_parcelname = row.names(regionalaverage.myelin))Myelin Maps
#Myelin Water Fraction Imaging
MWF.cifti <- read_cifti("/Volumes/Hera/Projects/corticalmyelin_development/Maps/MyelinWaterFraction_Liu2019/source-Liu2019_desc-MWFmap_space-fsaverage_den-164k_atlas-glasser360.pscalar.nii")
MWF <- data.frame(MWF = MWF.cifti$data, orig_parcelname = names(MWF.cifti$Parcel))
#Myelin Basic Protein Gene Expression
MBP.cifti <- read_cifti("/Volumes/Hera/Projects/corticalmyelin_development/Maps/AHBA_magicc/expression_maps/source-magicc_desc-MBPexpression_space-fsLR_den-32k.pscalar.nii")
MBP <- data.frame(MBP = MBP.cifti$data, orig_parcelname = names(MBP.cifti$Parcel))
MBP <- rbind(MBP, MBP) #reflect gene expression across hemispheres
MBP$orig_parcelname[1:180] <- glasser.regions$orig_parcelname[1:180]
#Sensorimotor-Association Axis
SAaxis.cifti <- read_cifti("/Volumes/Hera/Projects/corticalmyelin_development/Maps/S-A_ArchetypalAxis/FSLRVertex/SensorimotorAssociation_Axis_parcellated/SensorimotorAssociation.Axis.Glasser360.pscalar.nii")
SAaxis <- data.frame(SA.axis = rank(SAaxis.cifti$data), orig_parcelname = names(SAaxis.cifti$Parcel))df.list <- list(regionalaverage.myelin, MWF, MBP, SAaxis) #dfs to merge
myelin.maps <- Reduce(function(x,y) merge(x,y, all=TRUE, sort=F), df.list)
myelin.maps <- merge(myelin.maps, glasser.regions, by = c("orig_parcelname"), sort = F)
myelin.maps <- myelin.maps[!(myelin.maps$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),] #exclude low SNR parcelsR1 Indexes Cortical Myelin Content Across Cortical Regions
Whole-brain Map of R1
myelin.maps %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = mean.R1), , colour=I("gray50"), size=I(.06)) +
scale_fill_gradientn(colors = myelin.colorbar, limits = c(0.54, 0.58), oob = squish, na.value = "gray75") + theme(legend.position = "none")
R1 Correlation with Myelin Water Fraction Imaging
MWF Map
myelin.maps %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = MWF), colour=I("gray50"), size=I(.06)) +
scale_fill_gradientn(colors = myelin.colorbar, limits = c(0.02, 0.05), oob = squish, na.value = "gray75") + theme(legend.position = "none")
ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/MWF_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.5, height = 3)Correlation
##
## Spearman's rank correlation rho
##
## data: myelin.maps$mean.R1 and myelin.maps$MWF
## S = 2261576, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.6422757Correlation Plot
myelin.maps %>% filter(MWF < 0.08) %>% #remove 2 MWF major outliers in PHA1
ggplot(., aes(x = MWF, y = mean.R1)) +
geom_point(color = c("#8AABD5FF"), size = 0.2) +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
theme_classic() +
labs(x="\nMyelin Water Fraction", y="R1\n") +
scale_y_continuous(limits = c(0.525, 0.61), breaks = c(0.53, 0.55, 0.57, 0.59, 0.61)) +
theme(
axis.text = element_text(size = 6, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) 
R1 Correlation with Myelin Basic Protein Gene Expression
MBP Map
myelin.maps %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = MBP), colour=I("gray50"), size=I(.06)) +
scale_fill_gradientn(colors = myelin.colorbar, limits = c(-0.65, 1), oob = squish, na.value = "gray75") + theme(legend.position = "none")
ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/MBP_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.5, height = 3)Correlation
## Warning in cor.test.default(myelin.maps$mean.R1, myelin.maps$MBP, method =
## c("spearman")): Cannot compute exact p-value with ties##
## Spearman's rank correlation rho
##
## data: myelin.maps$mean.R1 and myelin.maps$MBP
## S = 3214088, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.4916123Correlation Plot
myelin.maps %>%
ggplot(., aes(x = MBP, y = mean.R1)) +
geom_point(color = c("#8AABD5FF"), size = 0.2) +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
theme_classic() +
labs(x="\nMyelin Basic Protein Expression (z)", y="R1\n") +
scale_y_continuous(limits = c(0.525, 0.61), breaks = c(0.53, 0.55, 0.57, 0.59, 0.61)) +
theme(
axis.text = element_text(size = 6, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) 
R1 Correlation with the Sensorimotor-Association Axis
S-A axis
myelin.maps %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = SA.axis), colour=I("gray50"), size=I(.06)) +
scale_fill_gradientn(colors = myelin.colorbar, trans = 'reverse', na.value = "gray75") + theme(legend.position = "none")
ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/SAaxis_corticalmap.pdf", device = "pdf", dpi = 500, width = 5.5, height = 3)Correlation
##
## Spearman's rank correlation rho
##
## data: myelin.maps$mean.R1 and myelin.maps$SA.axis
## S = 10280562, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## -0.6261257Correlation Plot
myelin.maps %>%
ggplot(., aes(x = SA.axis, y = mean.R1)) +
geom_point(color = c("#8AABD5FF"), size = 0.2) +
geom_smooth(method = "lm", color = c("black"), formula = y ~ x, linewidth = 0.2) +
theme_classic() +
labs(x="\nSensorimotor-Association Axis", y="R1\n") +
scale_y_continuous(limits = c(0.525, 0.61), breaks = c(0.53, 0.55, 0.57, 0.59, 0.61)) +
theme(
axis.text = element_text(size = 6, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) 
R1 Indexes Cortical Myelin Content Across Cortical Depths
#Compute mean R1 in each region at each cortical depth
regional_means <- function(input.df){
meanmap <- input.df %>% select(contains("ROI")) %>% colMeans() %>% as.data.frame() %>% set_names("mean.R1")
meanmap$orig_parcelname <- row.names(meanmap)
meanmap <- merge(meanmap, glasser.regions, by = c("orig_parcelname"), sort = F)
return(meanmap)
}
regional.myelin.alldepths <- lapply(myelin.glasser.7T, function(depth) {
regional_means(depth)
})Regional R1 Distributions by Depth
regional.myelin.alldepths <- do.call(rbind, regional.myelin.alldepths) #regional mean R1 at all depths
regional.myelin.alldepths$depth <- substr(row.names(regional.myelin.alldepths), 1, 7) #assign depth variable
regional.myelin.alldepths$depth <- factor(regional.myelin.alldepths$depth, ordered = T, levels = ordered_depths)
regional.myelin.alldepths <- regional.myelin.alldepths[!(regional.myelin.alldepths$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),] #exclude low SNR parcels
regional.myelin.alldepths %>%
ggplot(aes(x = mean.R1, y = depth)) +
geom_violin(aes(fill = depth), color = "white") +
stat_summary(fun = median, geom = "point", shape = 23, size = 2.25, color = "white") +
theme_classic() +
scale_fill_manual(values = depth_colorbar) +
scale_x_continuous(limits = c(0.495, 0.65), expand = c(0,0)) +
scale_y_discrete(expand = c(0, 0)) +
xlab("\nR1") +
ylab("Cortical Depth\n") +
theme(
legend.position = "none",
axis.text = element_text(size = 6, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) 
ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/Depthplot_regionalR1.pdf", device = "pdf", dpi = 500, width = 2.6, height = 4.9)for(this.depth in c(unique(regional.myelin.alldepths$depth))){
depth.plot <- regional.myelin.alldepths %>% filter(depth == this.depth) %>% filter(label != "lh_L_10pp") %>%
ggseg(.data = ., atlas = "glasser", mapping = aes(fill = mean.R1), colour=I("gray50"), size=I(.06)) +
paletteer::scale_fill_paletteer_c("grDevices::PuBuGn", direction = -1, limits = c(0.5, 0.62), oob = squish, na.value = "gray75") + theme(legend.position = "none")
print(depth.plot)
ggsave(filename = sprintf("/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/%s_corticalmap.pdf", this.depth), device = "pdf", dpi = 500, width = 4.75, height = 2.25)
}
R1 Myelin Profile in Individual Regions
Group Average Profiles
examplar.areas.alldepth <- rbind(area4.myelin.alldepths, area46.myelin.alldepths, areaa24.myelin.alldepths)
ggplot() +
geom_smooth(data = examplar.areas.alldepth, aes(x = mean.R1, y = depth, group = orig_parcelname, color= orig_parcelname), se = F, linewidth = .8) +
scale_color_manual(values = c(myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2], myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2])) +
theme_classic() +
xlab("\nR1") +
ylab("Cortical Depth\n") +
theme(
legend.position = "none",
axis.text = element_text(size = 6, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2))
ggsave(filename = "/Volumes/Hera/Projects/corticalmyelin_development/Figures/Figure1/Area_profiles_group.pdf", device = "pdf", dpi = 500, width = 2.4, height = 1.65)Individual Profiles
#R1 values in all regions at all depths for all participants
myelin.glasser.7T.alldepths <- do.call(rbind, myelin.glasser.7T)
myelin.glasser.7T.alldepths$depth <- substr(row.names(myelin.glasser.7T.alldepths), 1, 7) #assign depth variable
myelin.glasser.7T.alldepths$depth <- factor(myelin.glasser.7T.alldepths$depth, ordered = T, levels = ordered_depths)
#Wide to long format
cols_to_pivot <- names(myelin.glasser.7T.alldepths)[grepl("ROI", names(myelin.glasser.7T.alldepths))]
myelin.glasser.7T.alldepths <- myelin.glasser.7T.alldepths %>% pivot_longer(cols = all_of(cols_to_pivot), names_to = "orig_parcelname", values_to = "R1")
#Select exemplar regions of interest
myelin.glasser.7T.alldepths <- myelin.glasser.7T.alldepths %>% filter(orig_parcelname == "L_4_ROI" | orig_parcelname == "R_4_ROI" | orig_parcelname == "L_46_ROI" | orig_parcelname == "R_46_ROI" | orig_parcelname == "L_a24_ROI" | orig_parcelname == "R_a24_ROI")
#Format grouping variables
myelin.glasser.7T.alldepths$subject_id <- as.factor(myelin.glasser.7T.alldepths$subject_id)
myelin.glasser.7T.alldepths$orig_parcelname <- as.factor(myelin.glasser.7T.alldepths$orig_parcelname)
#Plot for the oldest 3 participants
myelin.glasser.7T.alldepths <- myelin.glasser.7T.alldepths %>% filter(age > 32)ggplot() +
geom_smooth(data = myelin.glasser.7T.alldepths, aes(x = R1, y = depth, color = orig_parcelname, group = interaction(orig_parcelname, subject_id)), se = F, linewidth = .4) +
scale_color_manual(values = c(myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2], myelin.colorbar[4], myelin.colorbar[3], myelin.colorbar[2])) +
theme_classic() +
xlab("\nR1") +
ylab("Cortical Depth\n") +
theme(
legend.position = "none",
axis.text = element_text(size = 6, family = "Arial", color = c("black")),
axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
axis.line = element_line(linewidth = .2),
axis.ticks = element_line(linewidth = .2)) 
Region-wise R1 Correlations Between Depths
regional.myelin.alldepths.wide <- regional.myelin.alldepths %>% pivot_wider(id_cols = orig_parcelname, names_from = depth, values_from = mean.R1) %>% select(-orig_parcelname)
regional.myelin.alldepths.wide <- regional.myelin.alldepths.wide %>% select(depth_7, depth_6, depth_5, depth_4, depth_3, depth_2, depth_1)
depth.R1.cors <- cor(regional.myelin.alldepths.wide) # compute correlation matrix
ggcorrplot(corr = depth.R1.cors, method = "square", type = "upper", show.diag = T, lab = F, lab_size = ) +
scale_fill_gradientn(colors = c("#edf6ff", "#c2d4ed", "#97adcc", "#345c91"))